#include "mark.h"
#include <cstdlib>
#include <cmath>
#include <fstream>
#include <ctime>

double reg::Rand(const double &min, const double &max) //Генерация равномерно распределенной случайной величины в диапозоне [min,max]
{
	return (double)rand()/(double)RAND_MAX*(max-min)+min;
}

void reg::variable(double *array)
{
	double x = 0, y = 0, s = 0;
	while((s > 1) || (s == 0))
		{
			x = Rand(-1,1);
			y = Rand(-1,1);
			s = x*x + y*y;
		}

	array[0] = x * sqrt(-2/s*log(s));
	array[1] = y * sqrt(-2/s*log(s));

}

void reg::variable(const double& M1,const double& Sig1,const double& cov1,
				   const double& M2,const double& Sig2,const double& cov2,double **array,const int& size) //Генерация нормально распределенной случайной величины
{

	double b[2];
	variable(b);
	array[0][0] = M1 + Sig1*b[0];
	array[1][0] = M2 + Sig2*b[1];
	for(int i = 1; i < size; i++)
	{
		variable(b);
		double ro1 = (cov1==0 ? 0 : cov1/(Sig1*Sig1)), ro2 = (cov2==0 ? 0 : cov2/(Sig2*Sig2));
		array[0][i] = M1 + array[0][i-1]*ro1 + b[0]*sqrt(Sig1*Sig1*(1-ro1*ro1));
		array[1][i] = M2 + array[1][i-1]*ro2 + b[1]*sqrt(Sig2*Sig2*(1-ro2*ro2));
	}
}

void reg::sample(const Matrix& X, double *array)
{
	double average = 0, square = 0, average_square = 0;
	double *pX = X.getpArray();

	for(int i=0; i < X.getSize(); i++)
		{
			average += pX[i];
			average_square += pX[i]*pX[i];
		}
	average /=X.getSize();
	average_square /=X.getSize();
	for(int i=0; i < X.getSize(); i++)
		square += (pX[i]-average)*(pX[i]-average);
	square /=X.getSize();
	array[0] = average;  //среднее значение выборки
	array[1] = square; //несмещенная дисперсия
	array[2] = sqrt(square); //стандартное отклонение
	array[3] = average_square; //среднее квадратов
}

void reg::mnk(const Matrix& Y, const Matrix& X, double *array)
{
	double average_xy = 0, a = 0, b = 0;
	double *pX = X.getpArray();
	double *pY = Y.getpArray();

	for(int i=0; i < X.getNumRow(); i++)
		{
			average_xy += pX[i]*pY[i]/X.getNumRow();
		}

	double c[4];//sample от x
	double d[4];//sample от y

	sample(X,c);
	sample(Y,d);

	a = (c[3]*d[0]-c[0]*average_xy)/(c[3]-c[0]*c[0]);
	b = (average_xy-c[0]*d[0])/(c[3]-c[0]*c[0]);

	array[0] = a;
	array[1] = b;
}

void reg::dmnk(const Matrix& C, const Matrix& I, const Matrix& G, const Matrix& Y, double *mark)
{
	// массив mark будет содержать коэффициенты c0, c1, i0, i1
	double r_CG[2];
	double r_IG[2];

	mnk(C,G,r_CG);
	mnk(I,G,r_IG);

	Matrix IkpG(G.getNumRow(),1);
	Matrix CkpG(G.getNumRow(),1);

	double *p_IkpG = IkpG.getpArray();
	double *p_CkpG = CkpG.getpArray();
	for(int i=0; i < G.getNumRow(); i++)
	{
		p_IkpG[i] = 0;
		p_CkpG[i] = 0;
	}

	double *pG = G.getpArray();

	for(int i=0; i < G.getNumRow(); i++)
	{
		p_IkpG[i] += pG[i] + (r_IG[0] + r_IG[1]*pG[i]);
		p_CkpG[i] += pG[i] + (r_CG[0] + r_CG[1]*pG[i]);
	}

	mnk(C,IkpG, r_CG);
	mnk(I,CkpG, r_IG);


	mark[0] = r_CG[0]; //a11
	mark[1] = r_CG[1]; //a12
	mark[2] = r_IG[0]; //a21
	mark[3] = r_IG[1]; //a22

	mark[1] /= (mark[1] + 1); //c1
	mark[0] *= (1 - mark[1]); //c0
	mark[3] /= (1 + mark[3]); //i1
	mark[2] *= (1 - mark[3]); //i0
}

void reg::MonteCarlo(Matrix& C, Matrix& I, const Matrix& G, Matrix& Y, const Matrix& Pi, const Matrix& B,
				const double& M1,const double& Sig1,const double& cov1,
				const double& M2,const double& Sig2,const double& cov2)//
{
	double *c = C.getpArray(),
		   *i = I.getpArray(),
		   *g = G.getpArray(),
		   *y = Y.getpArray();
	
	double **eps = new double*[2];
	eps[0] = new double[C.getNumRow()];
	eps[1] = new double[C.getNumRow()];

	variable(M1,Sig1,cov1,
				  M2,Sig2,cov2,eps,C.getNumRow());
				  
	Matrix X(2,1);
	Matrix PiX(2,1);
	double *x = X.getpArray();

	x[0] = 1;
	
	Matrix Binv = B;	
	Binv.inverse();
	
	Matrix Eps(2,1);
	double *veceps = Eps.getpArray();
	
	Matrix Igrek(2,1);
	double *igrek = Igrek.getpArray();
	for(int j=0; j < C.getNumRow(); j++)
	{
		x[1] = g[j];
		PiX = Pi*X;
		veceps[0] = eps[0][j];
		veceps[1] = eps[1][j];
		Igrek = PiX + Binv*Eps;
		c[j] = igrek[0];
		i[j] = igrek[1];
		y[j] = c[j] + i[j] + g[j];		
	}	
}

double reg::K(const double &u)
{
	return(1.0/sqrt(2*M_PI) * exp(-u*u*0.5));
}

double reg::MarkKernel(const double& sigma, const double& value, const double *const sample, const int& N)
{
double hN = 6.8*sigma*pow(N,-0.2);
double sum = 0;

for(int i = 0; i < N; i++)
	{
		sum += K((value - sample[i])/hN);
	}
return sum/N*hN;
}

double reg::min(const double* const sample, const int& N)
{
	double min=sample[0];
	for(int i = 1; i < N; i++)
	if(sample[i]<min) {min = sample[i];}
	return min;
}

double reg::max(const double* const sample, const int& N)
{
	double max=sample[0];
	for(int i = 1; i < N; i++)
	if(sample[i]>max) {max = sample[i];}
	return max;
}

void reg::printMarkKernel(const double* const sample, const int& size, const double& coldots, const std::string filename)
{
	double beg = reg::min(sample,size), end = reg::max(sample,size);
	double step = (end - beg)/coldots;
	std::ofstream out(filename.c_str());
	while(beg <= end)
	{
		out << beg << "\t" << reg::MarkKernel(sample[size],beg,sample,size) << std::endl;
		beg += step;
	}
	out.close();
}
